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The buckling of elastic bodies is a common phenomenon in the mechanics of solids. Wrinkling 
of membranes can often be interpreted as buckling under constraints that prohibit large amplitude 
deformation. We present a combination of analytic calculations, experiments, and simulations to 
understand wrinkling patterns generated in a bilayer membrane. The model membrane is composed 
of a flexible spherical shell that is under tension and that is circumscribed by a stiff, essentially 
incompressible strip with bending modulus B. When the tension is reduced sufficiently to a value 
0", the strip forms wrinkles with a uniform wavelength found theoretically and experimentally to be 
A = 2-k{B /aY^^ . Defects in this pattern appear for rapid changes in tension. Comparison between 
experiment and simulation further shows that, with larger reduction of tension, a second generation 
of wrinkles with longer wavelength appears only when B is sufficiently small. 

PACS numbers: 46.32.+X, 46.70.Hg, 46.25.-y 

I. INTRODUCTION 

Under excessive load a rigid body will buckle. A classic analysis of buckling dates to Euler who considered a 
thin rod under compressive stress along its axis. Above a critical stress the rod deforms with a wavelength twice 
as long as the rod, as the system relieves compressive energy at the expense of bending energy. Variations of Euler 
buckling appear in a range of geometries; however, under certain constraints, deviations from this behavior occur, and 
shorter wavelength patterns emerge. For example, rods subjected to a sudden compression buckle dynamically with 
a wavelength that scales inversely with the impacting stress [l| . A second example is the buckling of microtubules in 
biological cells for which the wavelength is reduced due to mechanical coupling to the surrounding cytoskeleton that 
energetically penalizes large amplitude deviations A related, well studied phenomenon occurs in the wrinkling 
of stiff elastic membranes attached to thick soft substrates d, 0, H, [E 0, Hi • Stretching of the substrate makes large 
amplitude, and hence long wavelength, deformations prohibitive, and the membrane wrinkles with a wavelength set by 
a compromise between the bending and stretching energies. A number of applications have been advanced that exploit 
the resulting periodic modulation such as in diffraction gratings Q , stretchable electronics [13, EH and metrological 
tools [T^ , indicating the potential of membrane wrinkling as an approach for patterning materials. 

In order to explore the variety of contexts in which one can controUably wrinkle membranes, we present a model 
of a bilayer membrane with a constraint that similarly prohibits long wavelength deformations. The membrane is 
composed of an inner layer that is a flexible spherical shell and an outer layer that is a stiff strip circumscribing the 
sphere. Initially, the sphere is stretched under an isotropic tension, while the strip layer is free of stress. When the 
tension on the sphere is reduced sufRciently, its radius decreases, and the strip buckles with a well defined wavelength. 
We have investigated this bilayer system in simple experiments by applying thin strips of paint around the equator of 
inflated rubber balloons, as shown in the photograph in Fig. 1(a). When the air pressure in the balloon is reduced, 
the strip wrinkles with a wavelength A that is much smaller than the radius of the balloon, as illustrated in Fig. 1(b). 
Comparing measurements of this wavelength with quantitative analytic results, we determine the dependence of A 
on other measurable parameters. Notably, we find that A depends on the stiffness of the strip but is independent 
of the elastic properties of the spherical membrane to lowest order, giving this geometry a potential advantage as a 
metrological approach for determining the bending modulus of thin stiff films. We further show how defects in the 
wrinkling pattern form when the pressure is reduced rapidly. We also introduce a numerical model for the system, 
which we employ to explore hierarchical structures that can emerge with large reductions in pressure. 

II. CALCULATION OF THE INSTABILITY WAVENUMBER 

In considering the energetics of the system, we approximate the sphere to be perfectly flexible, so that no bending 
energy is associated with its deformation, while assuming the strip is inextensible so that its length is conserved. We 
further parameterize the distance from the center of the sphere to the strip r in polar coordinates, r(0), and assume 
that the wrinkling instability forms a sinusoidal deformation in the strip, r{(j)) = f + acosn(f>, where n specifies the 
wavenumber. The objective of the calculation is to determine the wavenumber at the initial instability, where a 
grow infinitesimally from a = 0. In particular, we assume na <^ R for any finite n, where R is the initial radius of 



FIG. 1: (Color online) Experimental model of the bilayer membrane, which is comprised of a stiff strip of paint circumscribing 
an inflated rubber balloon, (a) Initially, the strip is free of stress and undeformed. (b) When the balloon is partially deflated, 
the strip wrinkles with a uniform wavelength. 



the balloon. In order to find the wavenumber, we need to calculate the total energy of the system to 0{a^). The 
inextensibility of the strip implies that 

2ttR = / dl= dtf) Vr^ + r'^ = 27rf + + Oia^). (1) 
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where dl = V r'^ + r'^ Acj) is the element of strip length, r' = dr/d(/), and we make use of the infinitesimal nature of a 
to expand the integrand. We thus obtain 

=i? + acosn0-— — + 0(a^). (2) 
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Upon reduction of the air pressure in the balloon, the relevant mechanical energy of the bilayer has three separable 
contributions. The first is the energy associated with bending the strip, which has the form 
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where B is the strip's bending modulus and k is the local curvature, 
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Inputting Eq. (|4]) in ([3]) and again expanding in powers of a, we obtain 

B f^"" , (r2 + 2r'2 - rr"f ttB irBa^ , , 
iJbending = ^ / d0^— = — + ^r^{n^ - \f + 0{a^). (5) 
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We note that the first harmonic [n — \) has zero bending energy because it amounts to a mere shift of the strip in 
space as long as a << i?, and the lowest harmonic affecting the shape of the strip is n = 2. 
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The second contribution to the energy arises from the inward force that the rubber at the edges of the strip appHes 
as the sphere deflates and its equihbrium radius decreases. This force per unit length (two-dimensional analog of 
pressure) is p = —2asma where a is the angle between the strip and the rubber membrane at the edge of the strip, 
as illustrated in the inset to Fig. 2(a), and cr is the tension of the balloon. The associated energy is 



-^pressure — pA — ^ 
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d(/) r((?!))^ = 27ri?^o-sinQ; - 7r(7i^ - l)a^cr sina + 0{a^), 
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where A is the area enclosed by the strip. Again, the first harmonic costs no energy because it does not affect the 
shape of the strip. We note that this expression is strictly correct only prior to the initial instability. After the 
instability p should acquire an dependence on that could affect the subsequent development of the instability but 
should not impact the instability wavenumber at the onset. 

If the total energy of the bilayer were only -Ebending + ^'prossuro, then the minimum energy solution for buckling would 
correspond to the lowest available harmonic, n — 2, and the instability would follow the analog of the Euler result 
for a circle: the strip would buckle into an ellipse. However, the elastic nature of the sphere changes this situation in 
a dramatic way because the additional stretching of the membrane introduced by a sinusoidal deformation penalizes 
large amplitude, long-wavelength deformations. To evaluate this elastic contribution, we consider the limit where the 
radius of the balloon is much larger than the width of the strip, so that the region of the rubber membrane affected 
by the strip can be treated essentially as a cylindrically symmetric system of radius f . The solution we will reach is 
consistent with this approximation. As long as the deformations of the membrane are small compared to the radius, 
the shape of the elastic membrane near the strip will obey Laplace's equation V^it = [l3|. In cylindrical coordinates 
this is 
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where z specifies the direction along the balloon perpendicular to the strip and u is the local deviation of the distance 
from the center to the surface of the membrane away from the average radius. The energy associated with the elastic 
deformation of the rubber is 
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where a is again the tension of the sphere. Using the ansatz u{z,(f>) — /(z)acosn0, as suggested by Eq. ^ (with 
/(O) = 1 and /(±cx)) — 0), Laplace's equation becomes 



dz2 



r/ = 0. 



(9) 



Hence f{z) = exp {—\nz\/f). That is, the deformation of the balloon decays exponentially as a function of distance 
from the paint strip with decay length f/n, and this result justifies the cylindrical approximation as long as the 
instability wavenumber n is sufficiently large Using this solution in the energy of Eq. ^ gives 



-^elastic = TTO-a^ |n| + 0{a^ 
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Therefore, the total energy we must consider is E{a,n) = £^bending + -Eprossurc + ^'clastic- Anticipating that the 
instability occurs at a high wavenumber we neglect the difference between n? — 1 and n^. In this approximation the 
total energy of the strip to ©(a^) is 



E{a,n) = (Bn'^/R^ - 2nVsina + 2a\n\) . 
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As a function of n, this energy has a minimum at n = no given by the condition dE/dn = 0. The instability occurs 
when E(a,no) — 0. Thus, we obtain the set of equations 



Bn^/B? - 2n^(7sina + 2an 
sin a + 2a 

Solving these, we find that the instability occurs at the wavenumber 
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The balloon pulls on the strip at an angle a such that 
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The angle is small if n 3> 1 . 



III. MEASUREMENTS OF INITIAL INSTABILITY 



To compare the wavenumber of the instability predicted by Eq. ()13p with that produced experimentally, we fabri- 
cated a set of bilayer membranes like those shown in Fig. 1. Nominally spherical balloons were inflated to an initial 
pressure of approximately 3100 Pa, corresponding to an initial radius of approximately 0.15 m. Paint (Krylon, "Fusion 
for Plastic" ) was applied to each balloon in a strip about the equator. The strip width w was kept sufficiently narrow, 
w = 5 mm, to neglect the curvature of the balloon across the width. To assure uniform application of the paint, each 
balloon was mounted on a turntable, and the paint was applied with an aerosol spray can positioned 0.3 m from the 
balloon, while the balloon rotated at an angular velocity of 60 RPM. Paper shields taped to the balloons exposed the 
5 mm wide region around the balloons' equator to the paint. To tune the thickness and hence the bending rigidity 
of the strip, the paint was applied in multiple thin coats, with a delay of approximately 2 minutes between each coat 
for the paint to set. The final thickness of the paint layers for different balloons was varied between 25 /im and 160 
/im, while the variability in thickness of any individual strip was less than 10%. 

After the painting, the balloons were left undisturbed for 18 hours to allow the paint to dry fully. To investigate 
the wrinkling instability, we measured the circumference of the strip to obtain R and then slowly released air from 
each balloon until deformations in the strip appeared like those in Fig. 1(b). The final air pressure Pj in the balloon 
was then measured. A measuring tape was carefully fastened to the balloon parallel to the wrinkled strip, and we 
determined the final balloon radius Rf, which was typically about 5% smaller than R. Digital photographs of the 
strips were taken, and image analysis was performed to extract the wavelength A of the undulations. Measurements of 
A were made at several points around the circumference of each balloon. From these results an average A as a function 
of thickness was obtained. The bending modulus of a paint strip can be related to its thickness h hy B = i2(i^u^) 
where E is the Young's Modulus and v is the Poisson ratio which through Eq. (fT^ leads to the prediction 
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where the tension is related to the final pressure by cr = RfPf/2. 

Figure 2 displays the results for A determined from a set of strips with different paint thicknesses. The wavelength 
is plotted in Fig. 2(a) against R^^h? /R^Pf on logarithmic scales to test the scaling relation predicted by Eq. (|15p . 
The error bars associated with the abscissa in Fig. 2 result primarily from the 10% variability in the paint thickness 
h, while the uncertainties in A reflect the spread in the measured values for each thickness, which again presumably 
results primarily from the variability in h for each balloon. As Fig. 2(a) demonstrates, the results for sufficiently large 

/^2,j3s 0.34±0.05 

h support the scaling prediction in Eq (I15p well, as evidenced by the result of a power-law fit, A \ ' 

shown by the solid line. Excluded from the fit are three data points at the smallest h, which deviate from this 
scaling. Observations of the onset of the instability indicated that the assumption of inextensibility of the strip was 
not valid at these smallest values of h. Specifically, we observed that the wrinkles annealed away moments after 
their formation, indicating that the thinnest paint layers are able to contract to relieve the stress. Relaxation of 
the inextensibility condition modifies the constraint on the modulation amplitude a, expressed in Eq. (2), ultimately 

changing the resulting scaling to A '-^ ( ) . No such contraction of the paint strip was visible at larger h 

where the results display good agreement with the predicted scaling. 

/Jj2 ^3^1/3 

In order to test the wavelength prediction in Eq. quantitatively, we display in Fig. 2(b) A versus ( -^xp- ) 
on linear scales, excluding the three points at smallest h. Assuming a Poisson ratio typical for an acrylic elastomer, 
V = 0.35, we extract from the fit the Young's modulus of the painted layer, E = 2.7 ± 0.3 x 10^ Pa, a very reasonable 
value for the paint. As a way to correlate this value with another measurement, we performed an additional experiment 
using the cantilever method to obtain the Young's modulus. The results of that experiment gave E = 2.9 ± 0.4 x 10® 
Pa, in good agreement with the value extracted from the wrinkling instability, further indicating the quantitative 
validity of Eq. 
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FIG. 2: (Color online) Results for the wrinkling patterns formed in paint strips of varying thickness following balloon deflation. 
In (a) the logarithm of the wavelength A is plotted against the logarithm of the quantity ^Tp^ • The line is the result of power 



law fit, A 



in good agreement with the predicted scaling, Eq. (|15|l . The data at the smallest h is excluded 



from the fit because the inextensibility condition is not obeyed for the thinnest paint layers. The inset shows a schematic cross 
section of the strip and adjoining rubber. The rubber on average pulls inward on the strip with a force per unit length p — —2a 

sina. In (b) A is plotted against I ^ p \ on linear scales. The line is the result of a linear fit from which the Young's 
modulus of the paint is extracted. 



We note that highly uniform wrinkling patterns, like those in Fig. 1(b), form reliably when the pressure in the 
balloon is decreased sufficiently slowly. In contrast, when the balloon is rapidly deflated, defects in the pattern 
appear. Figure 3(a) shows a photograph of the highly periodic modulation formed from slow deflation, while Fig. 3(b) 
shows a pattern formed by rapidly deflating a balloon at a rate dP/dt « 2000 Pa/s. In addition to regions with the 
wavelength expected from Eq. (I15|) . labeled Ai, the pattern displays regions with distinctly longer wavelength, such 
as A2. Similar defects appear commonly in wrinkling phenomena 0, IH. Their appearance in the strip following rapid 
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FIG. 3: (Color online) Upon rapid deflation, defects in the wrinkling pattern appear, (a) shows a photo of a uniform pattern 
resulting from slow deflation. The instability in (b) follows rapid deflation, which creates segments with equilibrium wavelengths 
Ai separated by regions of non-optimal wavelengths A2. 

deflation suggests that they are a consequence of non-simultaneity of wrinkle nuclcation in different regions around the 
circumference. When the incommensurate modulations in such different regions meet, the adjoining segments must 
adopt non-optimal wavelengths to match the regions smoothly. Unfortunately, difficulty in controlling quantitatively 
the deflation rate at high rates prevented us from investigating this effect in detail. 

IV. SECOND GENERATION WRINKLING: EXPERIMENTS AND MODELING 

As the pressure in the balloon is further reduced beyond the onset of the instability, the experiments display one 
of two effects on the wrinkling pattern. For large paint thickness and hence for large B, a uniform pattern persists, 
and the amplitude of the wrinkles grows until the pattern becomes highly corrugated, as shown in Fig. 4(a). For 
smaller B, a second larger-wavelength modulation emerges, as shown in Fig. 4(b). This second wavelength recalls the 
hierarchical, self-similar patterns observed recently in the wrinkling of stiff surfaces on soft thick substrates . 

To investigate the emergence of this second generation of wrinkles, we simulate the bilayer with a discrete model 
similar to one introduced previously to study membrane elasticity p^ . The membrane is approximated by two 
concentric circular chains of particles. Each chain has N particles, with ri^i and r2,i specifying the position of the i*'' 
particle in the outer and inner chains, respectively. The bending of the membrane is specified through unit vectors 
■hi for each particle on the outer chain, where the vector orientation is defined as normal to the line connecting the 
particle's two nearest neighbors on the chain. The membrane energy is defined as 

N 2 N 

^ = X! ■ + ~ ~ %")^ + X! ^3(1^1, i - r2.t\ - 03)^ (16) 

i i—1 j — 1 i—1 

where the first term sets the bending energy, with J determining the bending modulus, and the second term specifies 
the stretching energy, with spring constants Vi and V2 determining the compression modulus of the inner and outer 
chains, respectively. The third term couples each particle to its neighbor in the other chain through a spring force with 
constant V3, and the equilibrium bond length 03 sets the thickness of the membrane. The simulations begin with both 
chains forming undistorted circles, which minimizes both the bending and stretching energies, and with the inner circle 
having a radius Rq. The initial equilibrium bond lengths ai and 02 are set at values determined by Rq and 03 such 
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FIG. 4: (Color online) Photographs of paint strips on balloons deflated far beyond the onset of the wrinkling instability. For 
thin layers of paint and hence small bending modulus, the deformation remains uniform and becomes highly corrugated, as in 
(a). For thick layers and large modulus, a second, longer wavelength instability forms, as in (b). 

that the initial stretching energy is zero. To simulate the shrinkage of the inner layer, we decrease the equilibrium 
bond length in the inner chain in small steps 5a2 keeping all other variables constant. Typically, ^02/02 = 10~^, 
where is the initial bond length. After each step, the particle positions, ri^i and r2,i, are adjusted to minimize the 
energy, Eq. (fTB|) . using the conjugate gradient method [13 ■ After a certain number of steps, the circular symmetry 
is suddenly broken and the membrane deforms with a small amplitude, small wavelength pattern. This wrinkling is 
preferred energetically over deformation into an Euler-like ellipse because the small wavelength pattern enables the 
two chains to maintain a larger length difference, thus reducing their respective stretching energies, while avoiding 
a large separation between chains that would contribute a large energy from the term coupling the chains together. 
Thus, in the simulations the coupling between the inner and outer chains provides the constraint that prohibits large 
amplitude deformations. 

We find the wrinkling of the simulated membranes displays many of the salient features observed experimentally. 
For example, defects like those in Fig. 3(b) develop when the step size for the shrinkage 5a2/a^ is made large. The 
wavelength in the simulations is observed to scale with bending constant as A '~ J^/^, in contrast to the scaling 
predicted by Eq. ([15]) and observed in Fig. 2(a). We ascribe this difference to the extensibility of the outer chain in 
the simulations, which as described above is expected to alter the scaling. A detailed description of these findings 
will be presented elsewhere . Here, we focus on the behavior as the equilibrium radius of the inner chain is further 
reduced beyond the onset of the instability. In agreement with experiment, the amplitude of the wrinkles grows, 
and eventually a second generation of longer-wavelength wrinkles suddenly appears only if the bending constant 
is sufficiently small. The onsets of the first and second generation instabilities are easily identified by the sudden 
appearance of pronounced peaks in the Fourier power spectrum of the particle displacements. The inset to Fig. 5 
shows examples of such Fourier spectra shortly before and immediately after the formation of second generation 
wrinkles in a simulation. 

Figure [5] displays the total change in 02 normalized by at which the first and second generation wrinkle patterns 
appear in a series of simulations with varying bending constant J. Above a certain value of J, the second generation 
of wrinkles never appears, and instead the amplitude of the first generation continues to grow until the simulations 
break down. (Breakdown occurs when the inner chain passes through the outer chain, reminiscent of the experimental 
situation at large B in which the first generation modulation grows to the point that r((/>) is no longer single valued.) 
Strikingly, the strain amplitudes at which the first and second generation wrinkles appear display very similar quasi- 
power-law dependence on J. As suggested previously [SJ , the membrane with first generation wrinkles can be viewed 
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FIG. 5: Simulation results for the onset of the wrinkling instability. The relative strain Aa2/a2 at which the first generation 
(solid) and second generation (open) instabilities occur are displayed as a function of bending constant J. The simulation 
parameters employed were A'' = 200, as = 0.3, -Ro ~ 15, and Vj = 50,2,30 for j = 1,2,3 respectively. The inset shows the 
Fourier power spectrum of the particle displacements shortly before and immediately after the formation of second generation 
wrinkles in a simulation using the same parameters except = 600 and J — 0.5. The first generation power-spectrum peak 
(n « 27) has been scaled by 10^ for clarity. 



as a membrane with a larger effective thickness and larger effective bending modulus, which in turn wrinkles with a 
larger wavelength to create the hierarchical pattern. The similarity between the curves in Fig. 5 lends support to this 
picture, suggesting that the second generation instability can be mapped to the first through a simple rcnormalization 
of J that would lead naturally to self-similar structures. 



V. CONCLUDING REMARKS 



In summary, this simple system of a stiff strip on a spherical membrane has allowed us to investigate in detail 
a buckling instability in which a constraint enforces a small wavelength deformation. In this case, the constraint 
derived from the large energy cost of membrane stretching that would accompany an Euler-like, large amplitude 
buckling deformation. As a final perspective on the consequences of this constraint, we note briefly that one can 
formulate this contribution to the energy in terms of a long-range interaction between local deformations of the 
different parts of the strip. Specifically, above we analyzed the energy of small fluctuations of the local strip radius 
u{(j)) = r{(f)) — f for a given harmonic «((/>) = acosntf). In a more general case the total energy (|lip can be written in 
terms of the Fourier transform of u(l): 



/dk 
— {Bk^ /2 + k'^ - Ra sin a + a\k\) \uk 
27r 



(17) 



where / is the position measured along the strip. The first two terms - the bending and pressure energies - have 
regular Fourier transforms, whereas the elastic term has a kink singularity at fc = 0. Returning to real space we find 
that the bending and pressure energies become local terms containing spatial derivatives. The singular elastic term 
translates into a long-range interaction: 



E[u{l)] = / dl 



B fd^u\ (du 
_j _i?^sin« - 



dh / dh u{h)V{h~h)u{l2), 



(18) 
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where the interaction kernel is 

and the parameter t is an ultraviolet cutoff regularizing the short-distance behavior of the interaction. This iden- 
tification of the constraint with an effective long-range interaction suggests a method for mapping such wrinkling 
phenomena onto other condensed matter systems where competition between long-range and short-range interactions 
lead to pattern formation, such as for example in the creation of stripe domains in thin ferromagnetic films [T^. An 
interesting direction for further study would be to investigate whether such an approach can aid in understanding and 
predicting wrinkling behavior in less simple situations and further whether it might lend insight into cases in nature, 
such as the buckling of biological membranes during growth [Toj . in which the membranes are observed to wrinkle. 
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